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We discuss the problem of proton motion in Hydrogen bond materials with special focus on ice. 
We show that phenomenological models proposed in the past for the study of ice can be recast in 
terms of microscopic models in close relationship to the ones used to study the physics of Mott- 
Hubbard insulators. We discuss the physics of the paramagnetic phase of ice at 1/4 filling (neutral 
ice) and its mapping to a transverse field Ising model and also to a gauge theory in two and three 
dimensions. We show that Hg~0 and HO~ ions can be either in a confined or deconfined phase. We 
obtain the phase diagram of the problem as a function of temperature T and proton hopping energy 
t and find that there are two phases: an ordered insulating phase which results from an order-by- 
disorder mechanism induced by quantum fluctuations, and a disordered incoherent metallic phase (or 
plasma) . We also discuss the problem of decoherence in the proton motion introduced by the lattice 
vibrations (phonons) and its effect on the phase diagram. Finally, we suggest that the transition 
from ice Ih to ice XI observed experimentally in doped ice is the confining-deconfining transition of 
our phase diagram. 

PACS numbers: 05.30.-d, 77.80.-e. 71.30.+h 
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I. INTRODUCTION 

Hydrogen bonds (or H-bond) are ubiquitous in physics, 
chemistry, and biology. The nature of H-bonds was stud- 
ied in great detail by Pauling who predicted the mixed 
chemistry of H-bondsi, namely, H-bond shares charac- 
teristics of ionic and covalent bonds. On the one hand, 
a water molecule has a large dipole moment due to the 
electronegative character of the O atom and therefore 
water molecules are attracted to each other. This is the 
classical aspect of the H-bonding. On the other hand, 
the sigma bonding between H and O is strongly covalent 
with clear quantum mechanical nature. In the process 
of formation of an ice crystal electrons from the sigma 
bond can be shared by two water molecules leading to 
a strong link between them. Compton scattering experi- 
ments have confirmed the quantum nature of H-bonds in 
ice crystals^. It has been clear since the early experiments 
in H-bond systems that although the physics of electrons 
in ice is important the protons are actually responsible 
for the amazing electrical properties of ice—. The current 
understanding of the motion of protons in H-bond mate- 
rials is mainly based in a few phenomenological models. 
In this paper we discuss a qualitative microscopic model 
which captures the basic quantum mechanical correla- 
tions of the proton system. 

We show that the physics of protons in ice is a clear ex- 
ample of a strongly correlated problem very similar to the 
ones discussed in strongly correlated electron systems 5 
reminiscent of the physics of Mott insulators^. More 
specifically, the motion of protons in ice is hindered by 
strong constraints that forbid single proton hopping and 
only allows for collective ring-exchange- like motion. This 
sort of systems are known to be closely related to the 
physics of gauge theories. As we will see below, in the 



phase of ice in which the protons are ordered, in analogy 
with the problem of confinement of quarks in hadronic 
matter, pairs of defects (anions and cations) cost an en- 
ergy which is linear with the separation between them. 

Although there is a vast literature on physics of ice^, 
we focus only on the quantum motion of protons. The 
concepts introduced here can be easily extended to the 
study of many different problems in H-bond materials. 
Ice shows many different solid phases depending on tem- 
perature and pressure. These phases are essentially con- 
trolled by the geometry of the H-bonds that link different 
water molecules. Furthermore, a H atom is located asym- 
metrically relative to the two O atoms, forming a double 
well structure in the bond. Thus, a H atom can sit in 
any of the two sides a bond. Hence, we can think of ice 
as made out of protons, H + , moving in a crystal lattice 
made out of O -2 ions. Therefore, solids made out of H- 
bonds are expected to have strong electrical properties. 
For example, ice exhibits a high static permittivity com- 
parable with the one of liquid water, electrical mobility 
that is large when compared to most ionic conductors (in 
fact, the mobility is comparable to the electronic conduc- 
tion in metals). These are striking properties since in the 
solid phase only protons can move in an ice lattice as the 
electrons form a band insulator 3 . 

A striking feature of ice is its extensive classical en- 
tropy at low temperatures 4 . This large entropy implies 
a macroscopic number of classically degenerate states at 
low temperatures. This is a situation very similar to frus- 
trated magnetic systems, which the best representative 
is precisely called spin ice models. Our study, however, 
focus entirely on the proton motion in ice and not on 
magnetism. As we show, the proton motion in ice can 
be mapped into a problem with pseudo-spins, in close 
analogy to some frustrated magnets. 
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The most successful explanation for the physical be- 
havior of ice was given by the phcnomcnological work 
of Bernal and Fowled in 1933 that gave rise to the so- 
called Bernal- Fowler (BF) rules: (i) the orientation of 
H2O molecules is such that only one H atom lies be- 
tween each pair of O atoms; (ii) each O atom has two 
H atoms closer to it forming a water molecule. Rule (i) 
prevents situations in which a H-bond has two H atoms. 
Rule (ii) does not allow for configurations in which each 
O has more than two H close to it. As it was shown 
by Paulingi, using a classical counting argument, these 
rules can account for most of the entropy measured in 
the experiments, and predict an extensive entropy even 
at T = 0. A more detailed calculation taking into account 
ring exchange of protons in ice crystal brings this num- 
ber even closer to the measured valued. Pauling's sta- 
tistical model has been very successful in explaining the 
distribution of protons in ice and has been confirmed ex- 
perimentally in nuclear magnetic resonance (NMR) and 
neutron scattering experiments^. 

Although Pauling's calculation explains the arrange- 
ment of protons in ice it fails to describe its electrical 
properties. The reason for this failure is due to the fact 
that any motion of the protons under the BF rules re- 
quires an extremely correlated behavioB^. In order to 
explain the electrical behavior another phenomenologi- 
cal concept was introduced, namely, the concept of de- 
fects. Defects, by definition, are local violations of the 
BF rules. The defects associated with violation of rule 
(i) are called Bjerrum defects^ and the ones associated 
with violations of rule (ii) are called ionization defects^. 
With the concept of defects one is able to explain most 
of the electrical properties of icem. 

Protons can move by the rigid rotation of the water 
molecule, by thermal activation over an energy barrier 
from one side of the H-bond to another, or by quantum 
tunneling under the energy barrier between the two sides 
of the bond. At low temperatures, the rigid motion of the 
molecule and thermal activation are exponentially sup- 
pressed and only quantum tunneling is allowed. Quan- 
tum tunneling is possible because the wavefunction of 
the proton is extended from one side to another in the 
bond. Although the phenomenological theories of ice can 
account for a great part of the experimental data there 
are still many experiments that remain unexplained such 
as anomalies in the specific heat in purei^ and doped 
icoi&. Originally it was proposed by Onsager— that these 
anomalies could be explained by a ferroelectric transition. 
It turns out, however, that there is no evidence for any 
polar effect in ice. Quite to the contrary, the BF rules 
state that ice is a rather non-polar material. In this pa- 
per we show that by a quantum mechanical mechanism of 
order-by-disorder, an ordered phase of protons emerges 
at low temperatures. The melting of this ordered state 
can account for the specific heat anomalies. 

The classical ice model is defined on a pyrochlore lat- 
tice and its planar representation is the checkerboard 
lattice. The planar model is equivalent to the six ver- 



tex model what has been solved exactly by Liebi&. It is 
known from these studies that the phase diagram of the 
classical planar system has two phases, an ordered anti- 
ferroelectric phase and a line of critical points with exten- 
sive entropy at zero temperature. The anti-ferroelectric 
phase corresponds to a staggered arrangement of proton 
positions in the H-bonds. Its quantum version has been 
studied by Moessner, SondhiiS, and collaborators22i2I. 

In this paper, we consider the quantum version of the 
ice problem and show that quantum fluctuations stabi- 
lize the anti-ferroelectric phase in the absence of defects 
by a mechanism of order-by-disorder. To this end, we 
use a gauge theory description of the planar quantum ice 
problem which can be easily extended to higher dimen- 
sional systems. This gauge theory has also being used in 
similar models such as the quantum dimer models on the 
square lattice^. 

This paper is organized as follows: In the next section 
we discuss the characteristic energy scales of the ice prob- 
lem and propose the minimal model for proton motion in 
ice. Here we discuss in detail the analogy and connection 
to frustrated quantum magnets. Section IlIII contains the 
mapping of the planar ice model to a gauge theory. In 
Section ITVl we determine the ground state of the prob- 
lem in the neutral sector. Here we show how an ordered 
state of protons arises from the order-by-disorder mecha- 
nism and how the ionic defects are confined in the planar 
model. In Section E] we discuss the gauge problem in 
three dimensions and obtain the phase diagram includ- 
ing thermal effects. We show that the three dimensional 
case has a confining-deconfining transition even at T=0. 
In Section IVII we briefly discuss the effects of lattice vi- 
brations on the proton motion. We argue that the sec- 
ond order phase transitions obtained in the case without 
phonons can become first order. Section lYlIl contains our 
conclusions and the comparison between the theory and 
the experimental measurements in hexagonal (Ih) ice. 



II. THE MODEL 

Consider a lattice with protons living in the bonds and 
where each vertex is to be interpreted as an O atom as 
shown in Fig^ The lattices we are going o discuss in 
this papers are the planar square lattice, the cubic lattice 
and the pyrochlore lattice. The protons can occupy two 
positions in their respective link. Let us divide any of 
those lattices into sub-lattices 1 and 2, and for each link 
i define the sites il and i2 neighboring a vertex of sub- 
lattice 1 and 2, respectively. We can define the proton 
occupation number n.; ja , a = 1,2 for each proton site ia 
of the system. 

The main energy scales in this problem come from 
the Coulomb repulsion between the protons: the on-site 
Coulomb repulsion, Uh (the so-called Hubbard term), 
the Coulomb repulsion between protons in the same H- 
bond, Ub, and the Coulomb repulsion between protons 
around the same O vertex, Uo^. The physical situa- 
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tion of interest for ice corresponds to the case where 
Uh 3> Ub ^ Uo- Notice that this condition ensures 
that there is only one proton per lattice site, only one 
proton per bond, and for 1/4 filled ice it implies that 
there are only two protons around the O vertex. Thus, 
the electrostatic repulsion between the protons leads to 
the ice rules. 

Since we have considered interactions up to second 
nearest neighbors one might wonder whether we should 
not include interactions with even longer range. In fact, it 
is well-known that the water molecules interact through 
long-range dipole forces that in our picture would be 
represented by long-range interactions between H + ions. 
The same problem arises in the context of spin-ic o?' 2 ? . 
Monte carlo simulations^ for the frustrated lattices dis- 
cussed here have shown that the long-range dipolar in- 
teractions "average out" to zero at large distances. This 
result shows that the ice rules are still obeyed in the ab- 
sence of long-range interactions and that the degeneracy 
of the classical state is not lifted. Moreover, improved 
Monte Carlo simulations with "non-local" updates de- 
tected the appearance of long range order only at very 
low energies, much smaller than the energy scales dis- 
cussed here-^. Thus, from now on the effects of long- 
range interactions between second nearest neighbors will 
be disregarded. 

The Bjerrum and ionic defects mentioned earlier only 
occur if the protons hop from site to site. There are 
essentially only two types of hopping in this lattice: hop- 
ping on the H-bond with energy ts and hopping across 
the O vertex, to- Because the hopping energy is an ex- 
ponential function of the distance it is easy to show that 
ts ^> to- The presence of ts allows for the hopping of 
the protons around the O vertices leading to the creation 
of ionic defects. The hopping to-, on the other hand, 
allows for two protons in the same H-bond and there- 
fore can lead to Bjerrum defects. The values of ts and 
to are much smaller than in electronic systems because 
of the much larger proton mass. Thus, we expect that 
Uo 3> ts making the proton motion analogous to the 
motion of electrons in Mott insulators. 

In what follows we assume that Uh 3> Ub oo and 
tB 3> to —> so that only ionic defects are possible in 
the problem. It is also clear that in this limit the classical 
ground state obeys the ice rules and that magnetic phe- 
nomena associated with the proton spin does not play any 
role. Within the assumption that the on-site Coulomb re- 
pulsion between protons is very large, we restrain to the 
case where rn a — 0, 1, Vi and the proton spin degrees of 
freedom can be ignored. The Hamiltonian for the proton 
motion on neutral planar ice with just one proton per 
bond, that is, 



nu+rii2 = l, (2.1) 
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FIG. 1: Planar ice model: the vertices are O sites and the 
dots are H sites. The symbols are explained in the text. 

is written as, 

H = Up ^ (nuriji + n i2 n j2 ) 

(i.i> 

- t B Y J ( c li c ^ + h - c ) ■ ( 2 - 2 ) 

i 

We now make use of the constraint Eq. i|2.1[l , and define 
pseudo-spin operators associated to link i as (we use units 
such that Ub = 1 = K)'- 

Sf = - (4^2 + h. C )j , 

Si = \ ( c li c i2 - h-c)j , 

St = \ (cj lCil - c\ 2 c t2 ) , (2.3) 

which, due to condition Eq. Q2.1|) . obey the spin algebra 
[Sf, Sj] = iSijSf. In terms of these operators the Hamil- 
tonian of Eq. I|2.2|) is written as 

H = jJ2s;s;-r 1 £sf, (2.4) 

where J = 2Uo and T = 2ts- In the limit of T = the 
ground state is highly degenerate because there are many 
configurations of the pscudospins that give the same en- 
ergy. One possible configuration is shown in FigEl 

By a duality transformation, in which each link corre- 
sponds to a vertex of the square dual lattice, it is easy 
to show that Eq. l|2.4|) describes also the Ising model in 
a transverse field on a lattice where now the spins are 
defined on vertices. In the case of the planar ice system, 
this lattice correspond to the checkerboard lattice (see 
FigEJ), while for the pyrochlore and cubic lattices one 
obtains corner sharing tetrahedra and hexades, respec- 
tively. It is clear that the ground state of the Hamiltonian 
of Ea. H2.4fl with V — is an eigenstate of the operator 
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FIG. 2: A possible configuration of the pseudospins that sat- 
isfies the BF rules. 

S z . This state corresponds the proton localized every- 
where in either side of the link. In order to minimize the 
energy the distribution of protons is such that the only 
contributing configurations have the total magnetization 
of each vertex is equal to zero. For the pyrochlore and 
planar lattices, this means that each O has two protons 
close to it and two away. For the cubic lattice one would 
have three protons close to the O and three away. The 
number of configurations that satisfy the ice rules grows 
exponentially with the size of the system. Taking the pla- 
nar case as an example, there are (4/3) 37V ^ 2 of themi&. 



fects corresponds to the six vertex model (see FigQJa)). 
The full classical problem is a 14 vertex model that, as 
far as we know, has not been solved exactly^. 

Let us pick, for example, a particular configuration of 
the system that satisfies the ice rule and let us flip a single 
proton. This creates a pair of + (fi30 + ) and — (OH - ) 
ionic charges at two neighboring O sites, corresponding 
to the configurations FigQJc) and FigQJb), respectively. 
Strictly at Y = 0, one can separate the + and — defects 
along a zig-zag-like trajectory at no extra energy cost. 
There is however an entropic price for separating such 
defects that introduces an effective interaction between 
them. These charge defects play the role similar to holons 
in quantum dimer models^. Two natural questions arise 
in this context. The first one concerns possible lifting of 
the degeneracy of the ground-states by quantum fluctu- 
ations when ts oi T are finite. The second question is 
whether in the presence of these quantum fluctuations 
the ionic defects can be separately freely (unconfined) as 
in the classical case. 




(a) 



(b) 



(c) 



FIG. 3: Effective spin lattice for the ice problem in the Mott 
regime. The pseudospins interact with an antiferromagnetic 
exchange J. 

In the checkerboard language (see FigEJ a sublattice 
of plaquettes contains next nearest neighbor interactions 
("crossed plaquettes"). According to the BF rules the 
plaquettes have total S z magnetization equal to as 
shown in Fig^a). The configurations that contain de- 
fects violate the ice rules as shown in Figs^Jb) and (c). 
Classically, as it is well known, the system without de- 



FIG. 4: (a) Configurations of the pseudospins that obey the 
BF rules; (b) ionic defect with charge — (OH~); (c) ionic 
defect with charge + (H30 + ). 

To answer the first question, let us consider the op- 
erator P that projects into the subspace of states that 
satisfies the ice rules. Let PHP^ 1 denote the Hamilto- 
nian of Ea. (|2.4(l projected onto the ice rules sector, and 
let us treat the effects of finite Y in perturbation the- 
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ory. The effective Hamiltonian, for the planar and cubic 
lattices, to the lowest non- vanishing order is (S = a/2): 

□ 

where ij denote the links belonging to a crossed plaque- 
tte. To lowest order we find t ~ r 4 /J 3 . For the py- 
rochlore lattice, the lowest order effective Hamiltonian is: 

H ^Sf = - f E (^Jk^lm^mn^ni + h - c ) > ( 2 ' 6 ) 

ijklmn 

where now the interaction is around an hexagonal plaque- 
tte corresponding to the smallest loop in the pyrochlore 
lattice, and t oc T 6 /J 5 . The planar lattice with the 
Hamiltonian of Ea. (|2.5[) has been studied at zero tem- 
perature T = Oi^. It can be mapped onto a height 
model in which quantum fluctuations select the ordered 
flat state corresponding to a Neel order in the checker- 
board model. It is worth mentioning that the effective 
Hamiltonian of Ea. H2.5[) and Ea. (|2.tj[l also describe tun- 
neling within the low energy manifold of the XX Z spin 
1/2 model on the checkerboard. To see that, we start by 
writing the XXZ Hamiltonian: 

H=J2 Us?S* + ?f(SiSj + SrSf)) , (2.7) 

<ij> ^ ' 

where < ij > stand for all the couple of links belonging 
to the checkerboard lattice. We assume Jxy Jz- At 
zeroth order in J xy the ground state manifold is identical 
to the one mentioned before, with total S z magnetization 
per crossed plaquette equal to zero. Considering then the 
XY term in the first non- vanishing order in perturbation 
theory leads to Ea. (|2~5|) with t oc Jl y /J z and Ea. (|2~B|) 
with t oc J^y /J%, respectively. 



III. MAPPING OF THE PLANAR QUANTUM 
MODEL TO A LATTICE GAUGE THEORY 

In this section we map the ice problem onto a gauge 
theory2£iS&. Consider the O lattice which is shown 
schematically in Fig^ The H-bonds form links between 
the O atoms and to each link we can assign a value given 
by the pseudo spin defined in the previous section. Asso- 
ciated with these links we can assign states correspond- 
ing to the two configurations of the protons in each link 
which we will denote as | + 1/2) and | — 1/2). The ice 
rules imply 

ten 

where a z \ ± 1/2) = ±| ±1/2). This Ising gauge theory is 
equivalent to the one considered in Ref. 29], and as in 



that case, it has a local U(l) symmetry generated by the 
unitary operator 

Uv = JJ e to.Ei 6 .«f j (3.2) 

Civ 

with each a v associated to the vertex v being arbitrary. 
This approach has been applied with success in the case 
of the quantum dimer models. 

There is an alternative and complementary approach 
which works in an enhanced Hilbert space and that 
was applied originally to the quantum dimer model as 
well^S. Let us consider a square lattice is spanned by 
vectors r = (nei + me2)a where n and m are integers, 
ei = x, e2 = y and a is the lattice spacing. On each link 
(i.e., H-bond) we define a variable ^ M (r) with fi = 1,2 
which can be —1/2 or 1/2, which transfer like a vector: 
^_ M (r + e p ) = — £fj,(r). It will be convenient to define a 
new variable = + 1/2 which can take values of or 
1. Let us denote by A M the discrete derivative, 

A jU 4(r)=4(r + e jU )-4(r). (3.3) 

It is easy to see that the ice rules are equivalent to impose 
the condition that 

EA^(r) = 0. (3.4) 

Only configurations of the form Fig^a) are allowed for 
in. Therefore, Eq.(J22J reflects the BF rules. 

Eq. H3.4fl has the form of Gauss law in electrodynamics 
without external charges, V • E = 0. We will then assign 
to violations of the ice rules, such as those in Fig^c) 
the condition of 53«=i 2 A/i^ 1 ") = +1 f° r the forma- 
tion of a H^O and X)u=i 2^A( r ) = — 1 f° r the for- 
mation of HO~. Thus, if we allow violations of the ice 
rules by formation of ionization "defects" we must have 
^ ! 2 A M ^(r) = Q where Q plays the role of the "ef- 
fective" charge of the defect. In the quantum theory 
we define a Hilbert space of states |€ M (r)) which are the 
eigenstates of the operator E^r): 

£ M (r)|^(r))=^(r)|^(r)). (3.5) 

In the quantum theory Gauss law is a constraint in the 
space of state: 

E A M -E M (r)|Phys.) = 0, (3.6) 

which defines the physical states of the system. Let V 
be the operator canonically conjugated to E^, 

[e v (r'),E l ,(r)}=iS T yS^. (3.7) 

Thus, we see that E^ plays the role of the electric field 
while V plays the role of the vector potential in electro- 
dynamics. 
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Let us now define a ring exchange operator such that 
it maps configuration of electric fields satisfying the ice 
rule conditions. In this context this amounts to require 
that these operators, when acting on a plaquette, change 
the configuration in a manner consistent with Ea. H3.6f) . 
For each plaquette p we define the operator e 1 ®^ where 



(3.8) 



This operator is gauge invariant in the sense that it 
commutes with the generator of gauge transformations, 
Afj,Efj,(r). When acting on the links one gets, 



e ±l * (p) Mr)) = Mr)±l) 



(3.9) 



We are now going to relax the constraint that ^(r) 
can be or 1 by enlarging the Hilbert space to all integer 
values of ^(r). However, we will penalize energetically 
the values of £ M (r) that are not or 1 by adding an extra 
term to the Hamiltonian of the form: 



2g^ 



r./j 



(3.10) 



It is obvious that when g — > £^(r) can only be or 
1. Together with the constraint of Ea.fflty. Ecl iETTH 
defines the classical problem and the BF rules. Note 
that, for example, the Neel state corresponds to a con- 
figuration of electric fields in which the links (horizontal 
and vertical) form a staggered configuration of flux and 
1. These fluxes form stair like lines (right-up- right-up...) 
winding around the system. We also note that if the 1/2 
term in the kinetic energy Ea. H3.10|) where not present, 
this problem would be equivalent to compa ct elec trody- 
namics with matter fields studied in Refs. ,27 28] . The 
quantum kinetic energy of this problem is given by: 



(3.11) 



which is analogous to the magnetic energy in compact 
electrodynamics. 

In order to introduce ionic "defects" one has to intro- 
duce the charge Q into the problem. Associated with Q 
one defines the quantum operator n(r), which counts the 
number of defects at each vertex, and its conjugate <j>(r) 
such that 



(r'),n(r)]=iS r 



(3.12) 



where n(r) = 0, ±1, ±2, ... The Hamiltonian for the "mat- 
ter" field <j) contains a term for the energy required to 
create any one of these charges. For instance, one could 
write, 



(3.13) 



where Eq is the energy required to create a OH -Hjj~0 
pair. Observe that Ea. (|3.13|l tends to suppress defects. 



The motion of the defects is given by a kinetic energy 
term which is 

H c = cos (A„#r) - M (r)) , (3.14) 

where A ~ T is the coupling constant. This kinetic energy 
is gauge invariant since it commutes with the generator 
of local gauge transformations 



where 



U = e iS,a(r)G{r) ^ 



G(r)= Y, A^(r)-n(r) 



(3.15) 



(3.16) 



The gauge theory of the ice problem is described by 
the Hamiltonian H = Hq + Hj + Hq + He defined by 
Eqs. (|3~TU|) . l|3TTT)l . (|3~T3"J). and pi) . In the next section 
we are going to study the ground state of such a theory 
in the neutral secton^i. 



IV. GROUND-STATE SELECTION IN THE 
NEUTRAL SECTOR AND CONFINEMENT 

It is interesting to consider the analysis of the U(l) 
lattice gauge theory of quantum dimer models given in 
Refs. |al3Gl | in the context of our calculation. In the neu- 
tral sector, since the electric field has zero divergence, it 
can be written as: 



E lx v) = e llv [A„S(r) + B„(r)] 



(4.1) 



where S is an integer valued function defined on the dual 
lattice with periodic boundary conditions and B v is de- 
fined on the links of the dual lattice. Since the Gauss law 
tells us that 

only non-trivial topological configurations of B v have to 
be considered. Let us choose, in particular, B„ con- 
structed in the following way: let us define B^ for /j, = 2 
to take the alternating values and 1 on vertical links 
from row to row. For horizontal links with /j, = 1, will 
alternate from and —1 in such a way that each positive 
oriented arrow in the vertical direction meets a negative 
oriented arrow in the horizontal direction (see Fig[SJ|. By 
choosing such a B^ and choosing S = on the dual 
sites one gets for E^ one of the two Neel configuration 
(the other configuration being obtained by assigning a 
0, 1 staggered value to S). 

We can now write the path integral representation in 
discrete imaginary time in the form of the a 3-D discrete 
Gaussian model: 



J) 



2g 



[A,(S(r,j)-S )} 2 , (4.2) 
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□ n t 



FIG. 5: Background field B™ 



where j is the discrete time coordinate (e and A n are the 
"lattice spacing" and the discrete derivative in the imag- 
inary time direction, respectively) , and S° is a staggered 
configuration alternating between (in say, sub-lattice 
"A") and 1/2 (in sub-lattice "B"). While a non-zero t 
allows for fluctuations of the field S, only some kind of 
fluctuations can be made without paying the energy of 
the huge "surface tension" 1/g. More precisely, if at some 
point of sub-lattice "A" the value of S is 0, then the 
neighboring values of S (belonging to sub-lattice "B") 
can fluctuate between the values —1 and with no cost 
in surface tension. 

Imagine now that we consider the most general situ- 
ation for the electric field Ea. l|4.1|l in which the vector 
field has a non-trivial winding. The action for the 
discrete Gaussian model becomes : 



j) 



r.J 



+ ± ]T [A,(S(r,j)-S )-B,] 2 . (4.3) 

9 r,j,n=l,2 

Let us also assume that we choose the configuration for 
with the minimal number of links where it is non-zero. 
One can show that, for any configuration with non-zero 
B^, either there is no configuration of the S degrees of 
freedom that minimizes the surface tension term at each 
point, or at each point connected to a link where B^ 
is non-zero, there is a unique value for the field S(r,j) 
(modulo a global constant) that satisfies both periodic 
boundary conditions and a minimal surface tension. In 
the limit g — > the values of the field S associated to 
those points are not allowed to fluctuate. Then, the 
computation of the partition function of a discrete Gaus- 
sian model in a topological sector having non-zero B^ 
is penalized by the surface tension term or gives rise to 
the freezing of some of the plaquette degrees of freedom 
S(r, j) which gives a smaller contribution than the sector 
B^ = (since the configuration space of the former is a 
subset of the later). Keeping in mind these arguments, 



from now on, we are going to consider only the sector 
B, t = 0. 

To deal with the discrete variable S we make use of 
the Poisson summation formula: 



E /(»)= E / 

n—~00 771 — — oo 

and write the partition function: 



f{4>): 



(4.4) 



(4.5) 



where, 
S 



e 

2fJ 



\ i2irm(r,j)<j>(r,j) 



2tc 



£ [A„fo(r,j)-S°)] ! 



(4.6) 



p,= l,2 



We redefine: <p — <f> — S° , and write <fi as cf>A and (f>B for 
sub-lattices A and B, respectively. Following Ref. [32| . 
we work on the dilute gas approximation and keep only 
"monopoles" of charge ±1. In this case one finds an 
effective action of the form: 



S 



1 

2te 



A 



o(pB 



z cos(27r<^4) + z cos(27T0s)| , 



(4.7) 



where z is the monopole fugacity. In order to take the 
continuum limit, we write 



Xi + X2 



5B = Xl ~ X2 



(4.8) 



Note that a non-zero expectation value for \2 corre- 
sponds to the two Neel orders depending on the sign of 
< X2 > while < xi >= corresponds to a disordered 
state. Expanding in derivatives of the fields we get a 
continuum action of the form: 



S = \ <f 2 nh I \ [i^(Vxi) 2 + K T {dr X i) 2 + a xj] 



X sin(27rxi)sin(27TX2)} 



(4.9) 



Where V is the spatial two-dimensional gradient, K x , 
K T , a, and A are the parameters of the effective coarse 
grained model that in principle, but with a large amount 
of effort, can be calculated from the microscopic theory. 
Here the continuum limit has being taken with the pre- 
scriptions: 



ed T x(r,T) + 0(e 2 ), 
x(r,r)+ae M - V%(r, r) + 0(a 2 ) , 



A x(r,j) 

e — > J dr and a 2 — > J d 2 r 
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where a is the lattice spacing. We have not written the 
time and space variations of the field \2 because this field 
is massive and therefore the low energy physics is domi- 
nated by the x\ term. Because of that we can integrate 
out the \2 term in perturbation theory in A in order to 
generate an effective field theory for the field xi- 111 this 
case we obtain a relativistic-like field theory of the form: 



<S c ff 



• <!^;(v xi 



7C0s(4?rxi)} 



(drXl) 2 ] 



(4.10) 



where, K is the stiffness, and 7 = tt\ 2 / (2a) (we have set 
the velocity of propagation of the field equals to one, for 
simplicity). Eq. I|4.10|l describes a Sine-Gordon problem 
in 2 + 1 dimensions. Notice that in 2 + 1 dimensions the 
stiffness has dimensions of energy and, since the only en- 
ergy scale in this problem is the hopping t (since g — > 0), 
one concludes that K oc t. — Since we have assumed that 
7 <C 1 we can study this problem using a renormaliza- 
tion group argumen t 34 ' 35 , that is, we study the relevance 
of the cosine operator perturbatively by integrating high 
energy modes in a shell between A and A + dA where A 
is the ultraviolet cut-off of the theory. In doing that the 
coupling 7 renormalizes a.sM^: 



~d£ 



7 



(4.11) 



where d£ = dA/A. Hence, 7 is a marginally relevant 
coupling indicating that the field Xl is "frozen" at the 
minima of the potential in (|4.10|l . that is, at (xi) = ±1/4 
and a gap opens in the spectrum of Xi (thus, both xi and 
X2 are massive). The above result could be also derived 
directly from (|4.9() by minimizing the potential energy in 
order to find: 

(Xi) = Tl/4 ; <X 2 > = ±— cos(2^( X2 )). (4.12) 

a 

Thus, the energetic cost of instanton configurations in 
this systems forces it to freeze on one of the configurations 
that minimizes the potential. As mentioned before, the 
resulting ground-states present the Neel order. Thus, 
we have shown by a mechanism of order-by-disorder that 
the classical degeneracy is lifted by quantum fluctuations 
and select a state with Neel order. This conclusion is 
consistent with numerical simulations 20-21 . 

To understand what is the effect of a finite temper- 
ature, it is better to come back to the discrete action 
(|4.9() . We fist notice that the coarse grained system has 
effective Z2®U(\) symmetry, corresponding to changing 
X2 —> —X2 and xi ~~ > Xi + n /2 (where n is an integer). 
The ground state occurs with the simultaneous breaking 
of the two symmetries leading to a flat configuration of 
the fields. As the temperature is increased we expect the 
restoration of the Z2 symmetry and a rough phase with 
algebraic correlations for the xi field. The effective two- 



dimensional action at finite temperatures has the form: 

S = J d 2 rdrS^ [K(Vxi f + a(x2) 2 ] 

+ Asin(27rxi)sin(27rx 2 ) - 7 cos(4ttxi)} .(4.13) 

where /3 = 1/T is the inverse temperature. Clearly, 
the phase transition from a disordered high temperature 
phase to the ordered low temperature phase occurs with 
the breaking of the Z2 <8> U(l) symmetry. As usual, the 
critical temperature, T c , is proportional to the stiffness K 
and therefore, T c oc K oc i, consistent with the statement 
that t is the only physical energy scale in this problem. 

Equation (|4.1U|) also allows us to see explicitly how 
defects that are created over such ground state are con- 
fined by following Polyakov's argumentM. While this is a 
result that one could easily anticipated from the nature 
of the ground state itself, its description in terms of a 
gauge theory will be useful in higher dimensional models 
discussed in the next section. As we argued above, and 
as also considered in similar systems on Ref . [^2 > we l°°k 
for an effective action of the field xi which couples to the 
matter field (the defects). Integrating over the massive 
field X2 and expanding x = Xi~ < Xi > m Ea. (|4.10() wc 
obtain the effective action: 



S = J drdr 1 1 [(Vx) 2 + (d T x) 2 } 



+ m 2 x 2 } , (4.14) 



where m 2 rj 2it 2 j. These fluctuations are then massive 
and decay over a length scale ~ 1/m reflecting the fact 
that two-component Coulomb gas has Debye screening. 
This is also the confinement scale which defines the string 
tension a of the defects, that is, the confining potential 
has the form: V(R) = aR where a ~ m. This argu- 
ment follows from Polyakov's analysis^ of the standard 
compact QED in 2 + 1 dimensions which gives the same 
effective action as in Ea. (|4.14|) . At finite temperature, as 
long as the system remains in the fiat phase, the ionic de- 
fects will be confined. Increasing further the temperature 
will make the system enter the rough phase, at which the 
effective theory can be described by a massless theory for 
the field xi with U(l) symmetry^. Creation operators 
of ionic defects correspond to vertex operators of the dual 
fields (or dislocations in the surface roughening language) 
which now have algebraic decaying correlation functions. 



V. THREE-DIMENSIONAL MODELS 

In 3 + 1 dimensions the situation is more interesting 
since there is a confining-deconfining transition in the 
normal compact QED at T = 0. It is then in the cubic 
and pyrochlore lattice where we do have a chance to see 
deconfined ionic excitations even at T = 0. This argu- 
ment has been used recently by Hermele et al. in the 
context of the spin ice on the pyroclore lattice^,. The 
simplest three-dimensional lattice where we can apply 
this description is the cubic lattice. On this lattice, we 
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can in fact write the Hamiltonian with precisely the same 
form as in 2 + 1 dimensions. The main difference with 
the standard QED is again the term in Ea. H3.1()(l which 
penalizes configurations with £^^0,1. In Ref. [3t| . 
numerical evidence is provided that due to this term in 
Eq. (|3.10|) the defects are deconfined, which is to say that 
in the current model corresponds to a quantum disor- 
dered phase of ice (proton liquid). However, this analysis 
is based on a "naive" continuum limit of a lattice model 
and therefore is only applicable deep in the deconfined 
phase. It is known from early studies of compact QED 
in 3 + 1 dimensions that this theory has a deconfincment- 
confincment quantum phase transition at a critical value 
of the coupling constant t, the resonance amplitud e)?^?? . 
Qualitatively this phase transition is driven by the pro- 
liferation of "monopole loops". Notice that the theory 
described here, being perturbative in t, cannot describe 
the details of the quantum critical point that occurs at 
finite values of t. Nevertheless, these general arguments 
do predict the existence of a deconfinement-confinement 
quantum transition at finite t = t c and T = The ac- 
tual location of this quantum critical point may be con- 
trolled by additional operators not included here, such 
as the short range part of the dipole interaction, which 
affect the "electric" terms but leave the "magnetic" (flip) 
terms untouched. It is likely that the model we are con- 
sidering here, which has a single coupling constant t, with 
the electric terms playing the role of a constraint, may 
be in the deconfined (or "Maxwell" ) phase. 

The proton phases of ice are characterized by the en- 
ergy required to separate two defects with opposite sign. 
In the confined phase (the ordered phase of ice) the en- 
ergy grows linearly with the separation between the ions. 
Therefore, in this phase it is not possible to separate 
the H3O" 1 " and HO~ ions at arbitrary distance with the 
application of an external electric field. This state is a di- 
electric with a finite polarizability and it is an insulator. 
However, in the deconfined phase the effective interac- 
tion between the ions obeys V(R) oc 1/R and therefore 
decreases with the distance between the ions. Thus, in 
this case an applied electric field can separate the two 
H^"0 and HO~ ions at arbitrary distance leading to a 
conducting state. Hence, the phase transition is a metal- 
insulator transition 6 and therefore the conductivity of ice 
should behave quite different in the two different sides of 
the transition. We should note that in the metallic phase 
the conduction is due to the collective motion of protons. 
This is a correlated metal. 

In Fig we depict the phase diagram of ice as a func- 
tion of temperature and the proton hopping energy t. 
In Fig|Bfa) we show the phase diagram for planar ice 
(2 + 1 dimensions) for which we have shown that there 
is no deconfining transition as a function of t at T = 0. 
As we argued in the last section, there will be, however, 
a phase transition at finite temperature T c (T c cx t for 
t — > 0) between the disordered metallic phase and the 
ordered insulating phase. In the case at hand, the argu- 
ments presented above indicate the presence of a line of 



second order phase transition^. As we are going to see 
in the next Section, the coupling to phonons can change 
the transition to first order. In the case of the pyroclore 
or cubic lattice (3 + 1 dimensions) we have argued that 
there is a confining deconfining transition as a function of 
t, as shown in FigHJb). In this case, there is a quantum 
critical point (QCP) at some value of t (say, t c ). 



T c 



Deconfined Disordered 
Proton Plasma 



Confined Ordered 
Proton Insulator 



(a) 



T 



Deconfined Disordered 



Proton Plasma 



Confined Ordered 



Proton Insulator 



tc 

(b) 



FIG. 6: Phase diagram for protons in neutral ice. (a) Planar 
ice; (b) Three dimensional ice. 



VI. COUPLING TO LATTICE VIBRATIONS 

So far we have not discussed the problem of lattice vi- 
brations in ice. At finite temperatures one expects that 
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the thermal motion of the O atoms to affect the proper- 
ties of the H atoms. Vibrations must be important since 
the mass of the H atoms are just 16 times smaller than 
the O atoms and therefore the Born-Oppcnhcimer ap- 
proximation is not guaranteed. In order to incorporate 
the recoil of the O atoms due to the H motion we assume 
that the phonon coordinates are coupled to the local pro- 
ton density. In this case, one has to add an extra term 
to the Hamiltonian of the form: 



Hp — - K n i,a,<j<li,aM 




where <7i iQj0 are the local phonon coordinates, Pi iQj0 their 
canonical momentum, u>i M are the phonon frequencies, 
and M the ion mass. Eq. (j6.1|> can be reduced to the 
two level system in the same way Eq. lj2.2jl is reduced to 
Ea. H2.4|l . Using Eq. (j2.3jl and defining the relative and 
center of mass coordinates: 



•Ki.a — Qi,l,a Q_i.2,a i 
v _ Qi,l,a + Qi,2,a 



(6.2) 



the Hamiltonian Eq. (|6.1|) reduces to 

TT qz Pi,a , M ^i,o.xla 



n? g Mu ita xl a 

4M 4 



(6.3) 



in the basis of S z : Si\aj) 



<jj \crj). It is easy to see that 



the problem can be diagonalized by shifting operators: 













- b L - 








U 3,a 



(6.5) 



and the energy of the problem is given by: 
E[{a,n}] 



> U J 



{id) 



Tlia 



a 

Win 



(6.6) 



where the first term is just the energy of the classical 
state and nj >0 is the phonon number. So, the phonon 
frequencies remain the same but the atoms positions are 
shifted by Sxi ia oc KOi/ {Mujf a ). Since there is no pro- 
ton order when T = the average lattice shift is zero. 
However, for a small value of V we see from Fig. EI that 
<Ji acquires an average expectation value leading to an 
overall shift in the atom positions. At finite tempera- 
ture this shift will lead to a discontinuity in the specific 
heat as it is well known in the case of cooperative tran- 
sitions of this sort^. Thus, the phase diagram second 
order phase transition of Fig|S] can be modified substan- 
tially. In 2 + 1 dimensions the finite temperature phase 
transition can change completely to first order while in 
3+1 dimensions a tricritical point must appear at some 
temperature T* so that for T < T* the phase transition 
is second order and for T > T* the transition becomes 
first order. 



where pi is conjugated to Xi and is conjugated to Xi. 
Notice that the center of mass coordinate decouples from 
the proton motion in this case. 

It is interesting to rewrite the full Hamiltonian of the 
problem in terms of the pseudo-spin operators and stan- 
dard creation b\ and annihilation bi operator for the 
phonons &t] = Sij). Using Eq. lj2.4jl and proton- 
phonon coupling term Eo. (j6.3jl we find: 



(6.4) 



where Ai. Q = k/ WMcJ^a. Ea. (j6.4jl describes an Ising 
model in a transverse field coupled to a dissipative envi- 
ronment a la Caldeira-Leggett 40 . This model has been 
studied recently in the context of quantum phase tran- 
sitions in metallic magnetic systems^i*^ and has very 
interesting properties. 

Let us consider the classical case of T <c J when the ice 
rules are obeyed. In this case the problem can be solved 



VII. CONCLUSIONS AND EXPERIMENTAL 
REALIZATIONS 

In this paper we have considered the phase diagram of 
protons in ice. We have given a detailed picture of two- 
dimensional ice and provide a qualitative picture of the 
three-dimensional problem such as the pyrochlore lattice. 
We have studied the order-by-disorder effect and ground- 
state selection by quantum fluctuations, which in this 
case is the tunneling of protons between the two sites in 
a Hydrogen bond. 

In the absence of quantum tunneling, the number of 
low energy states grows exponentially with the size of 
the system and the entropy of the system is macroscopic. 
In this classical background, excited states correspond- 
ing to ionic defects can be created and separated with- 
out further cost in energy. As we showed, the presence 
of quantum fluctuations change considerably this sce- 
nario. The degeneracy of the system is lifted by the 
order-by-disorder mechanism leading to a well defined 
ground state. The technique used to understand the role 
of quantum fluctuations is based on a mapping of the 
lattice problem onto a gauge theory. 
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Within the gauge theory one can describe the behav- 
ior of ionic defects as confined or deconfined depending 
whether the ground state of the system is ordered or dis- 
ordered, respectively. In a deconfined state the ionic de- 
fects behave like as in a correlated metal, while in the 
confined phase the system behaves as in an insulator. 
For the two-dimensional model, quantum fluctuations se- 
lect a a two-fold degenerate antiferroelectric ground-state 
over which ionic excitations are confined. The situation 
is different in the three-dimensional case where topologi- 
cal excitations can lead to a deconfined disordered state 
even at T = 0. In such case, the interplay of quantum 
fluctuations and the tuning of the microscopic tunnel- 
ing parameter is expected to give rise to a confinement- 
deconfincment transition for ionic defects, with their cor- 
responding dramatic effects in the dielectric properties of 
the material. 

Furthermore, we have argued that at finite tempera- 
tures there will be a phase transition between ordered 
(confined) and disordered (deconfined) ice. The nature 
of this phase transition depends on the fugacity associ- 
ated with the topological defects (monopoles) and can 
be either first or second order. In the pyroclore lattice 
the most likely scenario is that there is a tricritical point 
separating a first order from a second order phase tran- 
sition (see FigOJb)) that should be observable. While 
similar problems have already being addressed in terms 
of spin models*t*2*, the ice scenario brings a new perspec- 
tive to this kind of phenomenology due the more acces- 
sible experimental possibilities than in frustrated spins 
systems— 

It is clear from Fig|(Jb) that the way to tune the phase 
transition is by changing the value of the proton hopping 
energy t. The value of t depends on the overlap of the pro- 
ton wavefunction between the two sites in the Hydrogen 
bond. Since the mass of the proton is large its wavefunc- 
tion is more localized than in the electronic case and we 
expect that t is exponentially sensitive on the changes in 
the distance between O atoms. One clear way to exper- 
imentally change t is then by applying hydrostatic pres- 
sure P to an ice crystal at fixed temperature T, starting 
from the disordered phase, and measuring the ice con- 
ductivity or dielectric response as a function of pressure. 
One expects that the proton ordering discussed here oc- 
curs at lower pressures than the well-known structural 
phase transitions of ice*. 

Another way to change the 0-0 distance is by using 
"chemical pressure"—. One can dope water with salts 
like KOH that while in solution become K + and HO~. 
When in the solid phase the K + ion becomes trapped 
into the "cages" of O atoms while the HO - ions are as- 
similated into the ice lattice. Thus, in this case, doping 
also introduces charges into the lattice structure and the 
ground state is no longer neutral. Furthermore, since the 
K+ ions spread randomly over the ice lattice, this kind 
of doping also introduces disorder in the system. The 
attractive interaction between the positive K + ions and 
the negative O -2 atoms of the ice structure lead to a lo- 



cal decrease of the volume of ice "cages" , to an average 
decrease in the H-bond distance, and to an average in- 
crease in t. So, in first approximation the introduction 
of KOH is somewhat equivalent to the horizontal axis in 
FigED However, because there is introduction of charges 
and disorder in the system, the horizontal axis in FigEI 
is only roughly the concentration of KOH. 

Doping experiments with KOH have been performed 
more than 20 years age** in order to understand the fa- 
mous 72 K anomaly observed in the temperature depen- 
dence of the specific heat of undoped ice— The specific 
heat of pure ice at low temperatures increases with tem- 
perature as T 3 which is the characteristic of phonons in 
the material. However, it has been known for a long time 
that hexagonal ice (ice-Ih) has a bump in the specific heat 
at T » 72 K that was a theoretical mysterj*— By doping 
pure ice with KOH it was shown by Kawada in the 1970's 
that the specific heat bump is actually a very slow phase 
transition into an ordered proton state that was called 
ice Xlaii*. Permittivity experiments in doped ice have 
confirmed this scenario— Specific heat measurements in 
KOH doped ice showed a strong and highly hysteretic 
first order phase transition with a substantial loss of en- 
tropy in the low temperature phase— Neutron scatter- 
ing experiments on single crystals of ice have confirmed 
the ordering transitional*. While the critical temper- 
ature is weakly dependent on the amount of KOH, the 
loss of entropy is dependent on the KOH concentration. 
Furthermore, in accordance with our discussion in the 
previous section, the first order phase transition is asso- 
ciated with lattice distortions as seen in recent neutron 
scattering experiments^*. 

Hence, our theory provides a possible theoretical expla- 
nation for the phase transition between ice-Ih and ice-XI. 
We believe that the first order phase transition observed 
in the experiments is the one described in Fig|Hfb) and it 
would be very interesting to find out whether the quan- 
tum critical point can be studied by further doping of 
KOH or by application of pressure. Our theory indicates 
that at T — > and t < t c the protons in ice-Ih would 
make an unique state of matter, namely, a quantum pro- 
ton liquid, with deconfined ion excitations. 
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